library('rjags')
library('coda')
library('jagshelper')
library('jagsUI')
library('mcmcplots')  # caterplot
library('dplyr')

maple_syrup_example

CONFIDENCE PROFILE EXAMPLE: Maple Syrup Urine Disease

Objective: To estimate the probability of retardation for a case of MSUD without screening (theta.n), & change in retardation rate associated with screening (e.d)

# Model

model_code <- "model
{

  #All data assumed to arise from binomial distributions with the appropriate parameters
   
   r.r ~ dbin(r, n.r)
   r.s ~ dbin(phi.s, n.s)
   r.n ~ dbin(phi.n, n.n)
   r.em ~ dbin(theta.em, n.em)
   r.lm ~ dbin(theta.lm, n.lm)  

  #Define functional relationships
    # Prob retardation for a case of MSUD who is screened
     theta.sm <- phi.s * theta.em + (1 - phi.s)*theta.lm

    # Prob retardation for a case of MSUD who is NOT screened               
     theta.nm <- phi.n * theta.em + (1 - phi.n) * theta.lm  

    # Expected retardation per 100000 newborns who are screened
     theta.s <- (theta.sm * r) * 100000

    # Expected retardation per 100000 newborns who are NOT screened                         
     theta.n <- (theta.nm * r) * 100000

   
  # Change in retardation rate associated with screening

  e.d <- theta.s - theta.n


  # Prior distributions - 'non-informative' Beta(1,1) priors

   r ~ dbeta( 1, 1)
   phi.s ~ dbeta( 1, 1)
   phi.n ~ dbeta( 1, 1)
   theta.em ~ dbeta( 1, 1)
   theta.lm ~ dbeta( 1, 1)

}
" %>% textConnection


# Data

data <- list( 
  n.r = 724262, 
  r.r = 7, 
  n.s = 276, 
  r.s = 253, 
  n.n = 18, 
  r.n = 8,  
  n.em=10, 
  r.em = 2, 
  n.lm = 10, 
  r.lm = 10)


# Starting/initial values



parameters_to_save <- c("e.d") # What else?


results <- jags(data = data,
            # inits = initial_values,
            parameters.to.save = parameters_to_save,
            model.file = model_code,
            n.chains = 1, # length(initial_values),
            n.adapt = 100,
            n.iter = 50000,
            n.burnin = 20000,
            n.thin = 2)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 5
##    Unobserved stochastic nodes: 5
##    Total graph size: 30
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                mean        sd       2.5%        25%        50%        75%
## e.d      -0.3411557 0.1665444 -0.7418316 -0.4318169 -0.3127903 -0.2207268
## deviance 20.1054097 3.2767312 15.6809410 17.6874476 19.4715124 21.8245152
##             97.5% Rhat n.eff overlap0         f
## e.d      -0.10076   NA     1        0 0.9999333
## deviance 28.26577   NA     1        0 1.0000000
plot(results)

hivepi_6

Simple version of the HIV Epi model with just 6 data points.

# Model

model_code <- "model{
#  SET PRIORS   
   a ~ dbeta( 1,2)            
   z ~ dbeta (1,1)      
   b <- z * (1-a)             #  sets constraint (1-a-b > 0)       
   c ~ dbeta (1,1)             
   d ~ dbeta (1,1)             
   e ~ dbeta (1,1)            

# VECTOR p[1:6] HOLDS THE EXPECTED PROBABILITIES FOR EACH DATA POINT
   p[1] <- a
   p[2] <- b
   p[3] <- c
   p[4] <- d
   p[5] <- (d*b + e*(1-a-b))/(1- a)
   p[6] <- c*a + d*b + e*(1-a-b) 
    
#  LIKELIHOOD AND DIAGNOSTICS
   for(i in 1: 6) {
       r[i] ~ dbin(p[i],n[i])                                                                                                      
       rhat[i] <- p[i] * n[i]                                                                                                                    
       dev[i] <- 2 * (r[i] * log(r[i]/rhat[i])  +  (n[i]-r[i]) * log((n[i]-r[i])/(n[i]-rhat[i])))    
    } 
    resdev <- sum(dev[])                                                                                                         
}
" %>% textConnection


# Data

data <- list( 
  r=c(11044, 12, 252, 10, 74, 254),
  n=c(104577, 882, 15428, 473, 136139, 102287)
)


# Starting/initial values: None specified

parameters_to_save <- c("a", "b", "c", "d", "e", "p", "r", "rhat", "dev", "resdev")


results <- jags(data = data,
            # inits = initial_values,
            parameters.to.save = parameters_to_save,
            model.file = model_code,
            n.chains = 1, # How many chains?
            n.adapt = 100,
            n.iter = 50000,
            n.burnin = 20000,
            n.thin = 2)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 6
##    Unobserved stochastic nodes: 5
##    Total graph size: 96
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                  mean           sd         2.5%          25%          50%
## a        1.057108e-01 9.478321e-04 1.038570e-01 1.050640e-01 1.057094e-01
## b        1.404450e-02 3.680695e-03 7.651316e-03 1.142074e-02 1.375050e-02
## c        1.715934e-02 8.759215e-04 1.547653e-02 1.656171e-02 1.714567e-02
## d        2.186721e-02 6.112196e-03 1.140236e-02 1.751401e-02 2.137004e-02
## e        2.456504e-04 1.226242e-04 2.264500e-05 1.533394e-04 2.460119e-04
## p[1]     1.057108e-01 9.478321e-04 1.038570e-01 1.050640e-01 1.057094e-01
## p[2]     1.404450e-02 3.680695e-03 7.651316e-03 1.142074e-02 1.375050e-02
## p[3]     1.715934e-02 8.759215e-04 1.547653e-02 1.656171e-02 1.714567e-02
## p[4]     2.186721e-02 6.112196e-03 1.140236e-02 1.751401e-02 2.137004e-02
## p[5]     5.799799e-04 6.292478e-05 4.629310e-04 5.360930e-04 5.776298e-04
## p[6]     2.332568e-03 9.657936e-05 2.146206e-03 2.265966e-03 2.331312e-03
## r[1]     1.104400e+04 0.000000e+00 1.104400e+04 1.104400e+04 1.104400e+04
## r[2]     1.200000e+01 0.000000e+00 1.200000e+01 1.200000e+01 1.200000e+01
## r[3]     2.520000e+02 0.000000e+00 2.520000e+02 2.520000e+02 2.520000e+02
## r[4]     1.000000e+01 0.000000e+00 1.000000e+01 1.000000e+01 1.000000e+01
## r[5]     7.400000e+01 0.000000e+00 7.400000e+01 7.400000e+01 7.400000e+01
## r[6]     2.540000e+02 0.000000e+00 2.540000e+02 2.540000e+02 2.540000e+02
## rhat[1]  1.105492e+04 9.912143e+01 1.086105e+04 1.098728e+04 1.105477e+04
## rhat[2]  1.238725e+01 3.246373e+00 6.748461e+00 1.007309e+01 1.212794e+01
## rhat[3]  2.647343e+02 1.351372e+01 2.387719e+02 2.555141e+02 2.645234e+02
## rhat[4]  1.034319e+01 2.891069e+00 5.393315e+00 8.284128e+00 1.010803e+01
## rhat[5]  7.895789e+01 8.566517e+00 6.302296e+01 7.298317e+01 7.863794e+01
## rhat[6]  2.385914e+02 9.878813e+00 2.195289e+02 2.317789e+02 2.384629e+02
## dev[1]   1.004898e+00 1.414663e+00 8.763011e-04 1.040125e-01 4.549952e-01
## dev[2]   8.655384e-01 1.224290e+00 1.000349e-03 8.649075e-02 3.943995e-01
## dev[3]   1.302654e+00 1.622315e+00 1.423556e-03 1.702855e-01 7.011324e-01
## dev[4]   8.329961e-01 1.178366e+00 9.743222e-04 8.366343e-02 3.804350e-01
## dev[5]   1.192679e+00 1.638148e+00 1.377183e-03 1.250219e-01 5.548485e-01
## dev[6]   1.413025e+00 1.419265e+00 4.871902e-03 3.434598e-01 9.987599e-01
## resdev   6.611791e+00 3.019223e+00 2.706447e+00 4.409680e+00 5.958335e+00
## deviance 4.697498e+01 3.019223e+00 4.306964e+01 4.477287e+01 4.632153e+01
##                   75%        97.5% Rhat n.eff overlap0 f
## a        1.063410e-01 1.075857e-01   NA     1        0 1
## b        1.634924e-02 2.214316e-02   NA     1        0 1
## c        1.773858e-02 1.888898e-02   NA     1        0 1
## d        2.574387e-02 3.510430e-02   NA     1        0 1
## e        3.350695e-04 4.810107e-04   NA     1        0 1
## p[1]     1.063410e-01 1.075857e-01   NA     1        0 1
## p[2]     1.634924e-02 2.214316e-02   NA     1        0 1
## p[3]     1.773858e-02 1.888898e-02   NA     1        0 1
## p[4]     2.574387e-02 3.510430e-02   NA     1        0 1
## p[5]     6.208112e-04 7.083781e-04   NA     1        0 1
## p[6]     2.396131e-03 2.526763e-03   NA     1        0 1
## r[1]     1.104400e+04 1.104400e+04   NA     1        0 1
## r[2]     1.200000e+01 1.200000e+01   NA     1        0 1
## r[3]     2.520000e+02 2.520000e+02   NA     1        0 1
## r[4]     1.000000e+01 1.000000e+01   NA     1        0 1
## r[5]     7.400000e+01 7.400000e+01   NA     1        0 1
## r[6]     2.540000e+02 2.540000e+02   NA     1        0 1
## rhat[1]  1.112083e+04 1.125098e+04   NA     1        0 1
## rhat[2]  1.442003e+01 1.953027e+01   NA     1        0 1
## rhat[3]  2.736708e+02 2.914192e+02   NA     1        0 1
## rhat[4]  1.217685e+01 1.660433e+01   NA     1        0 1
## rhat[5]  8.451662e+01 9.643789e+01   NA     1        0 1
## rhat[6]  2.450931e+02 2.584550e+02   NA     1        0 1
## dev[1]   1.334511e+00 5.004900e+00   NA     1        0 1
## dev[2]   1.143152e+00 4.395280e+00   NA     1        0 1
## dev[3]   1.830508e+00 5.695681e+00   NA     1        0 1
## dev[4]   1.106618e+00 4.239322e+00   NA     1        0 1
## dev[5]   1.603688e+00 5.808709e+00   NA     1        0 1
## dev[6]   2.070227e+00 5.161423e+00   NA     1        0 1
## resdev   8.155144e+00 1.424058e+01   NA     1        0 1
## deviance 4.851834e+01 5.460377e+01   NA     1        0 1
plot(results)

hivepi_6_xval

Simple version of the HIV Epi model with just 6 data points. (Yes, it has exactly the same comment as last example. But how is this one different?)

# Model

model_code <- "model
{
  #  SET PRIORS   
   a ~ dbeta( 1,2)            
   z ~ dbeta (1,1)      
   b <- z * (1-a)             #  sets constraint (1-a-b > 0)       
   c ~ dbeta (1,1)             
   d ~ dbeta (1,1)             
   e ~ dbeta (1,1)            

# VECTOR p[1:6] HOLDS THE EXPECTED PROBABILITIES FOR EACH DATA POINT
   p[1] <- a
   p[2] <- b
   p[3] <- c
   p[4] <- d
   p[5] <- (d*b + e*(1-a-b))/(1- a)
   p[6] <- c*a + d*b + e*(1-a-b) 
                                              


#  LIKELIHOOD AND DIAGNOSTICS
  for(i in 1:3) {       r[i] ~ dbin(p[i],n[i])   
       rhat[i] <- p[i] * n[i]               
       dev[i] <- 2 * (r[i] * log(r[i]/rhat[i])  +  (n[i]-r[i]) * log((n[i]-r[i])/(n[i]-rhat[i])))    
}                                                                                                   

 for(i in 5:6) {       r[i] ~ dbin(p[i],n[i])   
       rhat[i] <- p[i] * n[i]               
       dev[i] <- 2 * (r[i] * log(r[i]/rhat[i])  +  (n[i]-r[i]) * log((n[i]-r[i])/(n[i]-rhat[i])))    
}        
    dev[4]<-0
    resdev <- sum(dev[])        

#cross-validation of data point 4
    r.rep ~ dbin(p[4],n[4])
    p.xval <- step(r.rep - r[4])  -   0.5 * equals(r.rep,r[4])    
                                                                              
}
" %>% textConnection


# Data (same as before)

data <- list( 
  r=c(11044, 12, 252, 10, 74, 254),
  n=c(104577, 882, 15428, 473, 136139, 102287)
)


# Starting/initial values: None specified

parameters_to_save <- c("a", "b", "c", "d", "e", "p", "r", "rhat", "dev", "resdev", "p.xval")


results <- jags(data = data,
            # inits = initial_values,
            parameters.to.save = parameters_to_save,
            model.file = model_code,
            n.chains = 1, # How many chains?
            n.adapt = 100,
            n.iter = 50000,
            n.burnin = 20000,
            n.thin = 2)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 5
##    Unobserved stochastic nodes: 6
##    Total graph size: 93
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                  mean           sd         2.5%          25%          50%
## a        1.057004e-01 9.434816e-04 1.038805e-01 1.050636e-01 1.056868e-01
## b        1.358809e-02 3.974768e-03 6.948151e-03 1.070619e-02 1.324155e-02
## c        1.715228e-02 8.648453e-04 1.545648e-02 1.656562e-02 1.715344e-02
## d        2.158887e-02 1.514051e-02 9.007164e-04 9.683381e-03 1.975353e-02
## e        2.912632e-04 1.766201e-04 1.365644e-05 1.397866e-04 2.818106e-04
## p[1]     1.057004e-01 9.434816e-04 1.038805e-01 1.050636e-01 1.056868e-01
## p[2]     1.358809e-02 3.974768e-03 6.948151e-03 1.070619e-02 1.324155e-02
## p[3]     1.715228e-02 8.648453e-04 1.545648e-02 1.656562e-02 1.715344e-02
## p[4]     2.158887e-02 1.514051e-02 9.007164e-04 9.683381e-03 1.975353e-02
## p[5]     5.851106e-04 6.389430e-05 4.664383e-04 5.410095e-04 5.828977e-04
## p[6]     2.336223e-03 9.625111e-05 2.151411e-03 2.271167e-03 2.334889e-03
## r[1]     1.104400e+04 0.000000e+00 1.104400e+04 1.104400e+04 1.104400e+04
## r[2]     1.200000e+01 0.000000e+00 1.200000e+01 1.200000e+01 1.200000e+01
## r[3]     2.520000e+02 0.000000e+00 2.520000e+02 2.520000e+02 2.520000e+02
## r[4]     1.000000e+01 0.000000e+00 1.000000e+01 1.000000e+01 1.000000e+01
## r[5]     7.400000e+01 0.000000e+00 7.400000e+01 7.400000e+01 7.400000e+01
## r[6]     2.540000e+02 0.000000e+00 2.540000e+02 2.540000e+02 2.540000e+02
## rhat[1]  1.105384e+04 9.866647e+01 1.086351e+04 1.098723e+04 1.105241e+04
## rhat[2]  1.198469e+01 3.505746e+00 6.128269e+00 9.442856e+00 1.167905e+01
## rhat[3]  2.646254e+02 1.334283e+01 2.384626e+02 2.555744e+02 2.646433e+02
## rhat[5]  7.965637e+01 8.698506e+00 6.350044e+01 7.365249e+01 7.935512e+01
## rhat[6]  2.389653e+02 9.845237e+00 2.200614e+02 2.323109e+02 2.388288e+02
## dev[1]   9.934995e-01 1.418756e+00 8.343463e-04 9.890653e-02 4.442804e-01
## dev[2]   1.057669e+00 1.524609e+00 1.103423e-03 1.021089e-01 4.862977e-01
## dev[3]   1.277279e+00 1.540366e+00 1.826850e-03 1.727566e-01 7.207671e-01
## dev[4]   0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## dev[5]   1.298287e+00 1.741848e+00 1.261177e-03 1.391038e-01 6.313327e-01
## dev[6]   1.361105e+00 1.385083e+00 4.040936e-03 3.168897e-01 9.505058e-01
## resdev   5.987839e+00 2.894255e+00 2.402793e+00 3.877467e+00 5.322853e+00
## p.xval   4.479000e-01 4.855766e-01 0.000000e+00 0.000000e+00 0.000000e+00
## deviance 4.221527e+01 2.894255e+00 3.863022e+01 4.010490e+01 4.155028e+01
##                   75%        97.5% Rhat n.eff overlap0 f
## a        1.063193e-01 1.075926e-01   NA     1        0 1
## b        1.601359e-02 2.237791e-02   NA     1        0 1
## c        1.773533e-02 1.885851e-02   NA     1        0 1
## d        3.036678e-02 5.680448e-02   NA     1        0 1
## e        4.358797e-04 6.170680e-04   NA     1        0 1
## p[1]     1.063193e-01 1.075926e-01   NA     1        0 1
## p[2]     1.601359e-02 2.237791e-02   NA     1        0 1
## p[3]     1.773533e-02 1.885851e-02   NA     1        0 1
## p[4]     3.036678e-02 5.680448e-02   NA     1        0 1
## p[5]     6.267139e-04 7.150763e-04   NA     1        0 1
## p[6]     2.399721e-03 2.528530e-03   NA     1        0 1
## r[1]     1.104400e+04 1.104400e+04   NA     1        0 1
## r[2]     1.200000e+01 1.200000e+01   NA     1        0 1
## r[3]     2.520000e+02 2.520000e+02   NA     1        0 1
## r[4]     1.000000e+01 1.000000e+01   NA     1        0 1
## r[5]     7.400000e+01 7.400000e+01   NA     1        0 1
## r[6]     2.540000e+02 2.540000e+02   NA     1        0 1
## rhat[1]  1.111855e+04 1.125171e+04   NA     1        0 1
## rhat[2]  1.412399e+01 1.973732e+01   NA     1        0 1
## rhat[3]  2.736206e+02 2.909490e+02   NA     1        0 1
## rhat[5]  8.532020e+01 9.734978e+01   NA     1        0 1
## rhat[6]  2.454602e+02 2.586357e+02   NA     1        0 1
## dev[1]   1.310196e+00 5.060578e+00   NA     1        0 1
## dev[2]   1.392216e+00 5.282334e+00   NA     1        0 1
## dev[3]   1.830402e+00 5.563814e+00   NA     1        0 1
## dev[4]   0.000000e+00 0.000000e+00   NA     1        1 1
## dev[5]   1.751668e+00 6.159568e+00   NA     1        0 1
## dev[6]   1.969325e+00 4.995338e+00   NA     1        0 1
## resdev   7.366281e+00 1.337691e+01   NA     1        0 1
## p.xval   1.000000e+00 1.000000e+00   NA     1        1 1
## deviance 4.359371e+01 4.960434e+01   NA     1        0 1
plot(results)